エクセルVBAで高速フーリエ変換

Sub FFT()

' 2012.9.27 移植
' FFT(高速フーリエ変換)

Dim cr(100, 10000) As Variant, ci(100, 10000) As Variant
Dim co(10000) As Variant, si(10000) As Variant
Dim pi As Variant
Dim n(10000) As Variant
Dim l As Variant 'ポイント数(2の乗数)
Dim m As Variant
Dim fs As Variant 'サンプリング周波数[Hz]

pi = Application.WorksheetFunction.pi()

Call data_input(cr, ci, fs, l, m) 'データ入力
Call bit_reverse1(co, si, pi, l, n) 'nのビットリバースとωの作成
Call butterfly(cr, ci, co, si, l, m, n) 'バタフライ演算
Call bit_reverse2(cr, ci, l, m, n) 'スペクトルのビットリバース
Call spectrum(cr, ci, fs, pi, l) 'スペクトルの出力

End Sub

'--------------------- データ入力 --------------------------
Sub data_input(cr, ci, fs, l, m)
Do While Cells(i + 2, 2) <> ""
cr(0, i) = Cells(i + 2, 2).Value
i = i + 1
Loop
l = i
m = l
fs = CDbl(l - 1) / CDbl(Cells(l + 1, 1) - Cells(2, 1))
Cells(2, 11) = l
Cells(2, 10) = fs
i = 0
Do While 2 ^ i < l
i = i + 1
Loop
l0 = l
l = 2 ^ i
m = l
For i = l0 + 1 To l
cr(0, i) = 0
Next i
Cells(2, 12) = l

End Sub
'----------------- nのビットリバースとωの作成 ---------------
Sub bit_reverse1(co, si, pi, l, n)

For k = 0 To l - 1
' Let w(k) = complex(Cos(2 * (k) * pi() / l), -Sin(2 * (k) * pi() / l))
co(k) = Cos(2 * k * pi / l)
si(k) = Sin(2 * k * pi / l)
For j = 0 To log2(l) - 1
Let n(k) = n(k) + (Int(k / 2 ^ j) Mod 2) * 2 ^ (log2(l) - 1 - j)
Next j
Next k

End Sub
'---------------------- バタフライ演算 ----------------------
Sub butterfly(cr, ci, co, si, l, m, n)

For j = 0 To log2(l) - 1
Let kk = 2 ^ j - 1
Let m = m / 2
For jj = 0 To kk
For k = 0 To m - 1
cr(j + 1, k + 2 * jj * m) = cr(j, k + 2 * jj * m) + cr(j, k + 2 * jj * m + m) * co(n(2 * jj)) + ci(j, k + 2 * jj * m + m) * si(n(2 * jj))
ci(j + 1, k + 2 * jj * m) = ci(j, k + 2 * jj * m) - cr(j, k + 2 * jj * m + m) * si(n(2 * jj)) + ci(j, k + 2 * jj * m + m) * co(n(2 * jj))
cr(j + 1, k + 2 * jj * m + m) = cr(j, k + 2 * jj * m) + cr(j, k + 2 * jj * m + m) * co(n(2 * jj + 1)) + ci(j, k + 2 * jj * m + m) * si(n(2 * jj + 1))
ci(j + 1, k + 2 * jj * m + m) = ci(j, k + 2 * jj * m) - cr(j, k + 2 * jj * m + m) * si(n(2 * jj + 1)) + ci(j, k + 2 * jj * m + m) * co(n(2 * jj + 1))
Next k
Next jj
Next j

End Sub
'----------------- スペクトルのビットリバース ----------------
Sub bit_reverse2(cr, ci, l, m, n)
For k = 0 To l - 1
cr(log2(l) + 1, n(k)) = cr(log2(l), k)
ci(log2(l) + 1, n(k)) = ci(log2(l), k)
Next k
End Sub
'--------------------- スペクトルの出力 ----------------------
Sub spectrum(cr, ci, fs, pi, l)
For k = 0 To l - 1
Cells(k + 2, 4).Value = k * fs / l '周波数
Cells(k + 2, 5).Value = (2 / l) * cr(log2(l) + 1, k) '実部
Cells(k + 2, 6).Value = (2 / l) * ci(log2(l) + 1, k) '虚部
Cells(k + 2, 7).Value = (2 / l) * Sqr(cr(log2(l) + 1, k) ^ 2 + ci(log2(l) + 1, k) ^ 2) '振幅スペクトル
Cells(k + 2, 8).Value = (180 / pi) * atan2(cr(log2(l) + 1, k), ci(log2(l) + 1, k)) '位相スペクトル
Next k

End Sub
Function log2(x As Variant) As Variant
log2 = Application.WorksheetFunction.Log(x, 2)
End Function
Function atan2(x As Variant, y As Variant) As Variant
atan2 = Application.WorksheetFunction.atan2(x, y)
End Function